setwd("~/Probably/Downloads")

library(texreg)
library(car)
library(ggplot2)
library(texreg)
options(scipen=120,digits=2)
library(extrafont)
font_import()


data <- read.csv("Epperly_JLC2018_Aydin_replication_data_with_unincluded_type2_states.csv")

names(data)
# The file contains the data used by Aydin, as well as the additional following variables:
# miscoding, which is my flag for whether it was coded as a democracy but is only a type2 "democracy" according to the ACLP data.
# also included is included_by_aydin, which flags whether the observation was included in Aydin's paper (as I include those typ2 regimes that for some reason are not included in the original analysis)


# Fixing missing data in replication file (doing it here rather than in the file for transparency purposes)
data$system[data$countries=="Serbia"] <- 1
data$system[data$countries=="Taiwan"] <- 0
# Taiwan has no ppp pc gdp because data uses WB data. 
# CIA world factbook:
data$gdpcapitappp[data$countries=="Taiwan"] <- (17400+17200+18000+23400+25300+27500+29500+30100+31100)/9
data$loggdp[data$countries=="Taiwan"] <- log((17400+17200+18000+23400+25300+27500+29500+30100+31100)/9)


data2 <- subset(data,miscoding==0)
dim(data2)
# data2 is those countries data analyzes that are NOT type2 regimes, but in factt democracies
data3 <- subset(data,miscoding==1)
dim(data3)
# data3 is those countries included by data but that are actuallty type2 regimes
data4 <- subset(data,included_by_aydin==1)
dim(data4)
# data4 is the subset of data including those analyzed by data (minus the one missing country in the replication file provided by data)
data$miscoding_as_factor <- as.factor(data$miscoding)
data$miscoding2 <- recode(data$miscoding_as_factor,"NA='3';else=0")
data5 <- subset(data, miscoding==1 | miscoding2==3)
dim(data5)
# data5 is all type2 regimes, including those used in Aydin's article and those she omitted
data6 <- subset(data, included_by_aydin==0)
# data6 is the type2 regimes not included by data


plot1 <- ggplot(data2, aes(x=demqualityfinal)) + 
	stat_density(geom="line", size=1)
plot1 <- plot1 + stat_density(data=data3, geom="line", linetype="dashed",size=1)
plot1 <- plot1 + stat_density(data=data6, geom="line", linetype="dotted",size=1)
plot1 <- plot1 + scale_x_continuous(name='',limits=c(0,1))
plot1 <- plot1 + scale_y_continuous('Density',limits=c(0,4))
plot1 <- plot1 + theme_bw()
plot1 <- plot1 + theme(plot.title=element_text(size=25))
plot1 <- plot1 + theme(axis.title.y = element_text(size=16,angle=90,family="Garamond"),
          axis.text.y=element_text(size=16,family="Garamond"))
plot1 <- plot1 + theme(axis.title.x = element_text(size=16,vjust=0,family="Garamond"),
          axis.text.x  = element_text(size=16,family="Garamond"))
plot1 <- plot1 + annotate("text", x=0.82, y=3.65, label="Democracies",family="Garamond",size=8)
plot1 <- plot1 + annotate("text", x=0.82, y=3.4, label="(as in Cheibub et al.)",family="Garamond",size=6)
title1 <- substitute(paste(bold("A  "),"Level of democracy"))
plot1 <- plot1 + ggtitle(title1) 
plot1 <- plot1 + theme(title = element_text(size=18,family="Garamond"))
plot1 <- plot1 + theme(axis.title.x = element_text(size=14,family="Garamond"),
	axis.text.x=element_text(size=14,family="Garamond"))
plot1 <- plot1 + annotate("segment", x=0.55, xend=0.65, y= 3.65, yend= 3.65,size=1)
plot1

plot2 <- ggplot(data2, aes(x=polcomp)) + 
	stat_density(geom="line", size=1)
plot2 <- plot2 + stat_density(data=data3, geom="line", linetype="dashed",size=1)
plot2 <- plot2 + stat_density(data=data6, geom="line", linetype="dotted",size=1)
plot2 <- plot2 + scale_x_continuous(name='',limits=c(0,1))
plot2 <- plot2 + scale_y_continuous('Density',limits=c(0,4.5))
plot2 <- plot2 + theme_bw()
plot2 <- plot2 + theme(plot.title=element_text(size=25))
plot2 <- plot2 + theme(axis.title.y = element_text(size=16,angle=90,family="Garamond"),
          axis.text.y=element_text(size=16,family="Garamond"))
plot2 <- plot2 + theme(axis.title.x = element_text(size=16,vjust=0,family="Garamond"),
          axis.text.x  = element_text(size=16,family="Garamond"))
plot2 <- plot2 + annotate("text", x=0.25, y=4.25,
	 label="type2",family="Garamond",fontface="italic",size=8)
plot2 <- plot2 + annotate("text", x=0.51, y=4.25,
	 label="regimes (included)",family="Garamond",fontface="plain",size=8)
plot2 <- plot2 + annotate("segment", x=0.02, xend=0.18, y= 4.25, yend= 4.25,
	size=1, linetype="dashed")
plot2 <- plot2 + annotate("text", x=0.25, y=3.75,
	 label="type2",family="Garamond",fontface="italic",size=8)
plot2 <- plot2 + annotate("text", x=0.553, y=3.75,
	 label="regimes (not included)",family="Garamond",fontface="plain",size=8)
plot2 <- plot2 + annotate("segment", x=0.02, xend=0.18, y= 3.75, yend= 3.75,
	size=1, linetype="dotted")	
title2 <- substitute(paste(bold("B  "),"Political competition"))
plot2 <- plot2 + ggtitle(title2) 
plot2 <- plot2 + theme(axis.title.x = element_text(size=14,family="Garamond"),
	axis.text.x=element_text(size=14,family="Garamond"))
plot2 <- plot2 + theme(title = element_text(size=18,family="Garamond"))
plot2


res1 <- lm(judicial.independence.de.facto ~ polcomp * demqualityfinal,data=data4)
res2 <- lm(judicial.independence.de.facto ~ polcomp * demqualityfinal + common_law + 				as.factor(system) + juddejure + log(gdpcapitappp),data=data4)
res3 <- lm(judicial.independence.de.facto ~ polcomp * demqualityfinal,data=data2)
res4 <- lm(judicial.independence.de.facto ~ polcomp * demqualityfinal + common_law +
		as.factor(system) + juddejure + log(gdpcapitappp),data=data2)

screenreg(list(res1, res2, res3,res4),dcolumn=TRUE,
		custom.coef.names=c("Intercept","Political competition","Democracy level",
		"Competition*Democracy level","Legal system",
		"Semi-presidential system","Parliamentary system",
		"De jure independence","GDP per capita (ppp, logged)"),
		caption="Replication of Aydin (2013)",
		caption.above=TRUE,label="table:data",booktabs=TRUE
		)
# THIS ABOVE GETS YOU TABLE 2 in my article (Table 1 is just a table listing countries, not done in R)

# Code below produces the marginal effects plots for Figure 2 (a) and (b)

library(MASS)
beta <- coef(res2)
covvar <- vcov(res2)
# quartile 1 is approx. equal to Honduras, Kenya, Nicaragua
# quartile 3 is approx. equal to US, New Zealand, France
polconrange <- seq(min(data4$polcomp),max(data4$polcomp),length.out=100)
x <- cbind(1,polconrange,
	summary(data4$demqualityfinal)[2],1,0,0,1,
	mean(log(data4$gdpcapitappp),na.rm=TRUE),
	polconrange*summary(data4$demqualityfinal)[2])
set.seed(1)
beta.sim <- mvrnorm(10000,beta,covvar)
xb <- x %*% t(beta.sim)
confinf <- apply(xb,1,quantile,probs = c(0.025, 0.5, 0.975))
dat <- data.frame(x = polconrange, lower = confinf[1,], 
	middle = confinf[2,], upper = confinf[3,])
dat
x2 <- cbind(1,polconrange,
	summary(data4$demqualityfinal)[5],1,0,0,1,
	mean(log(data4$gdpcapitappp),na.rm=TRUE),
	polconrange*summary(data4$demqualityfinal)[5])
xb2 <- x2 %*% t(beta.sim)
confinf2 <- apply(xb2,1,quantile,probs = c(0.025, 0.5, 0.975))
dat2 <- data.frame(x = polconrange, lower = confinf2[1,], 
	middle = confinf2[2,], upper = confinf2[3,])

plot3 <- ggplot(dat, aes(x = x))
plot3 <- plot3 + geom_ribbon(aes(ymin = lower, ymax= upper),
		alpha=0.5,fill="gray70")
plot3 <- plot3 + geom_ribbon(data=dat2,aes(ymin = lower, ymax= upper),
		alpha=0.5,fill="gray70")
plot3 <- plot3 + geom_line(data=dat2,aes(y = middle), size=1, colour="black",linetype = 2)
plot3 <- plot3 + geom_line(aes(y = middle), linetype = 1,size=1)
plot3 <- plot3 + scale_x_continuous('Political competition',limits=c(0,1)) 
plot3 <- plot3 + scale_y_continuous('E(Judicial independence)',limits=c(2.5,8.75))
plot3 <- plot3 + theme_bw()
plot3 <- plot3 + theme(axis.title.y = 
	element_text(size=14,angle=90,family="Garamond"),
	axis.text.y=element_text(size=14,family="Garamond"))
plot3 <- plot3 + theme(axis.title.x = element_text(size=14,family="Garamond"),
	axis.text.x=element_text(size=14,family="Garamond"))
title3 <- substitute(paste(bold("A  "),"Including ",paste(italic("type2 ")),"regimes"))
plot3 <- plot3 + theme(axis.title.x = element_text(size=14,family="Garamond"),
	axis.text.x=element_text(size=14,family="Garamond"))
plot3 <- plot3 + ggtitle(title3) 
plot3 <- plot3 + theme(title = element_text(size=18,family="Garamond"))
plot3
# Produces Figure 2(a)

polconrange2 <- seq(min(data2$polcomp),max(data2$polcomp),length.out=100)
beta2 <- coef(res4)
covvar2 <- vcov(res4)
# quartile 1 is approx. equal to Ukraine, Armenia, or El Salvador
# quartile 3 is approx. equal to Denmark, Portugal, of Germany (and 0.03 higher than the US)
x3 <- cbind(1,polconrange2,
	summary(data2$demqualityfinal)[2],1,0,0,1,
	mean(log(data2$gdpcapitappp),na.rm=TRUE),
	polconrange2*summary(data2$demqualityfinal)[2])
beta.sim3 <- mvrnorm(10000,beta2,covvar2)
xb3 <- x3 %*% t(beta.sim3)
confinf3 <- apply(xb3,1,quantile,probs = c(0.025, 0.5, 0.975))
dat3 <- data.frame(x = polconrange2, lower = confinf[1,], 
	middle = confinf[2,], upper = confinf[3,])
x4 <- cbind(1,polconrange2,
	summary(data2$demqualityfinal)[5],1,0,0,1,
	mean(log(data2$gdpcapitappp),na.rm=TRUE),
	polconrange2*summary(data2$demqualityfinal)[5])
xb4 <- x4 %*% t(beta.sim3)
confinf4 <- apply(xb4,1,quantile,probs = c(0.025, 0.5, 0.975))
dat4 <- data.frame(x = polconrange2, lower = confinf2[1,], 
	middle = confinf2[2,], upper = confinf2[3,])

plot4 <- ggplot(dat3, aes(x = x))
plot4 <- plot4 + geom_ribbon(aes(ymin = lower, ymax= upper),
		alpha=0.5,fill="gray70")
plot4 <- plot4 + geom_ribbon(data=dat4,aes(ymin = lower, ymax= upper),
		alpha=0.5,fill="gray70")
plot4 <- plot4 + geom_line(data=dat2,aes(y = middle), size=1, colour="black",linetype = 2)
plot4 <- plot4 + geom_line(aes(y = middle), linetype = 1,size=1)
plot4 <- plot4 + scale_x_continuous('Political competition',limits=c(0,1)) 
plot4 <- plot4 + scale_y_continuous('E(Judicial independence)',limits=c(2.5,8.75))
plot4 <- plot4 + theme_bw()
plot4 <- plot4 + theme(axis.title.y = 
	element_text(size=14,angle=90,family="Garamond"),
	axis.text.y=element_text(size=14,family="Garamond"))
plot4 <- plot4 + theme(axis.title.x = element_text(size=14,family="Garamond"),
	axis.text.x=element_text(size=14,family="Garamond"))
title3 <- substitute(paste(bold("B  "),"Excluding ",paste(italic("type2 ")),"regimes"))
plot4 <- plot4 + theme(axis.title.x = element_text(size=14,family="Garamond"),
	axis.text.x=element_text(size=14,family="Garamond"))
plot4 <- plot4 + ggtitle(title3) 
plot4 <- plot4 + theme(title = element_text(size=18,family="Garamond"))
plot4
# Produces Figure 2(b)

